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Abstract 

Various properties of the general two-center two-electron integral over the explicitly correlated 
exponential function are analyzed for the potential use in high precision calculations for diatomic 
molecules. A compact one dimensional integral representation is found, which is suited for the 
numerical evaluation. Together with recurrence relations, it makes possible the calculation of the 
two-center two-electron integral with arbitrary powers of electron distances. Alternative approach 
via the Taylor series in the internuclear distance is also investigated. Although numerically slower, 
it can be used in cases when recurrences lose stability. Separate analysis is devoted to molecu- 
lar integrals with integer powers of interelectronic distances r\i and the vanishing corresponding 
nonlinear parameter. Several methods of their evaluation are proposed. 
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I. INTRODUCTION 



In order to achieve high accuracy for the nonrelativistic energy as well as for relativis- 
tic and quantum electrodynamics (QED) corrections in molecular systems, well optimized 
basis functions have to be applied. The explicitly correlated Gaussian (ECG) functions do 
not satisfy the cusp condition, therefore their use in the evaluation of higher order QED 
corrections [l| is problematic. Nevertheless, the recent calculations of the leading 0(a 3 ) 
QED effects performed for H 2 [2| using ECG functions, lead to the most accurate to date 
theoretical predictions of about cm -1 uncertainty for dissociation energies. To obtain 
even more accurate results and to include higher order QED effects we aim to use the basis 
of explicitly correlated exponential functions of the form 

<P — e ' 12 ' IA ' IB r 2A r 2B- \ l ) 

The notation for interparticle distances is explained after Eq. fl4]). The use of such exponen- 
tial functions in molecular applications is quite limited. Starting from the pioneering work 
of Kolos and Roothan [3j, the Neumann expansion of r^ 2 (at oj = 0) in spherical oblate 
coordinates |4( has been applied most often, see the most recent review by Harris in Ref. 
and the collection of works in Ref. [6|. This expansion and the more general one of 



e -u>ri2 r -i j n terms of spheroidal functions |4|, has been applied by several authors |7H11] in 
their accurate calculations for the hydrogen molecule. 

An alternative approach to perform integrals with exponential function in Eq. ([!]) has 



originally been proposed in Ref. [12[. Authors have obtained an analytic, although quite 
complex formula for the general four-body integral with exponential functions. By taking 
the inverse Laplace transform [see Eq. (^} in one of the nonlinear parameters, one can in 
principle obtain the general two-center two-electron integral in Eq. (T5]). Due to the very 
complicated analytic structure of the four-body formula, this inverse Laplace transform 
has not been applied so far. In our recent work [ljj, we have reformulated the problem 
of the calculation of the inverse Laplace transform, into the solution of some differential 
equations. From these differential equations satisfied by the master integral f(r) (see Eq. (T2J) 
below), one derived analytic recursion relations for integrals with positive integer powers of 
interparticle distances, assuming that the nonlinear parameter corresponding to ri2 vanishes, 
i.e. for standard molecular integrals. Using these recursions one obtained analytic results for 
integrals with James- Coolidge and extended Heitler-London basis functions. The application 



of analytic formulas has been demonstrated by the calculation of Born- Oppenheimer energies 



for H2 in Ref. [14( with accuracy of 10 -15 au, and for HeH + in Ref. 15| with accuracy of 
about 10~ 12 au. While James-Coolidge and Heitler-London basis functions work very well 
for the ground states, they are not equally good for excited states, where arbitrary values 
of nonlinear parameters are needed. This is the main subject of this work, to develop a 
computational technique for general two-center two-electron integrals, which can be used 
for the calculation of nonrelativistic energies, and also for relativistic and QED effects in 
two-electron diatomic molecule. Finally we think that various properties of the general two- 
center two-body integral derived here, can be applied for the calculation of nonrelativistic 
energies of an arbitrary diatomic molecule. This subject has recently been pursued also by 



Lesiuk and Moszyhski in Ref. 



16]. 



The general structure of this work is the following. In Sec. II we derive the differential 
equation for the master integral f(r) in Eq. fl2]). In Sec. Ill we find a solution in terms 
of one- dimensional integral over elementary functions and additionally consider few special 
cases. In Sec. IV the Taylor series of f(r) in the internuclear distance r is derived, and 
this shows the general analytic properties of the master integral /. In Sec. V integrals 
with positive powers of interparticle distances are obtained by differentiation of the master 
integral with respect to the corresponding nonlinear parameter. In Sec. VI we consider 
standard molecular integrals, i.e., with the vanishing nonlinear parameter corresponding to 
ri2. In Sec. VII we present a summary with possible further applications. 



II. THE MASTER INTEGRAL 



We introduce here notation and definitions following the previous work in Ref. [13|, and 
obtain the fourth order differential equation which is satisfied by the master integral f{r). 
Matrix elements of the nonrelativistic Hamiltonian between functions of the type ([I]) can 
be expressed in terms of /(r) and its derivatives with respect to nonlinear parameters. The 
master two-electron two-center integral f(r) is defined by 

or, / Or, e e e e e 

4tt J 4tt r 12 r 1A r 1B r 2A r 2B 

and the related class of integrals with the nonnegative integer n is 

/"d 3 ri f d 3 r 2 1 e -^r 1A e -u 2 r 1B e -w 2 r 2A e -w 3 r 2B 

/(r , n) = r y_y ______ — _ (3) 



where 1, 2 are indices of the electrons, A, B that of the nuclei, and r = tab is the distance 
between the nuclei. The notation for nonlinear parameters comes from the general four-body 
integral g defined by 

/ x [ d 3 Ol f d 3 2 f d 3 p 3 e - W lP^- W 2P2-mP3-U lP 23-U2P31-U 3 pi2 

9[Ui) = / - — / - — / -: , 4 

J 47T J 4tt J 4vr P23 P31 Pl2 Pi P2 P3 

with pi = fi2, f>2 = r 2 A, f>3 = r 2 B- Function g is related to / by a Laplace transform, namely 



d 3 r „, s e 



t r 



9(t)= / i-/(r)^-= / drf(r)e- tr } (5) 



o 



4 7r J w r 2 

and the opposite relation is the inverse Laplace transform 



i oo+e 

tr 



f(r) = — / dte tr g(t). (6) 

ATll J _i OQ+fL 



It was shown in Ref. 



131 ]. that the function g satisfies the first order differential equation in 



any of it's parameter a = Ui, wf. 

do 1 da „ 

or equivalently 

Vj^(Vjg) + P a = 0, (8) 

with the inhomogeneous term P a presented in Appendix A, and a being the sixth order 
polynomial in six variables 

222, 222, 222, 222,22/2, 2 2 2 2 2\ 

(J = u 1 u 2 w 3 + U 2 u 3 W 1 + U l U 3 W 2 + W 1 W 2 W 3 + U l W 1 [U 1 + W 1 — U 2 — u 3 — W 2 — w 3 ) 

,2 2 / 2,2 2 2 2 2\,22/2,2 2 2 2 2\ /n\ 

+U 2 U7 2 (« a + W 2 - W x - U 3 - W l - «J 3 ) + U 3 W 3 (U 3 + W 3 - U 2 - U x - W 1 ~ W 2 ). (9) 

For example, Eq. ([Zj) for a — w\ reads 

<y ^ h o a — 9 + P(wi, Ui, w 2 , u 2 ; w 3 ,u 3 ) = , (10) 

ow\ 2 aw i 

where the inhomogeneous term P(w\,u\\ w 2 ,u 2 ;w 3 ,u 3 ) is given by Eq. flA2p . The solution 
of this differential equation is presented in the ingenious work of Fromm and Hill [12] by 
direct integration of Eq. ([2]) in the momentum representation. A more compact formula was 



obtained by Harris 



17|. We do not present here their results, as its explicit form is quite 



long and we will not need it. 
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To obtain the fundamental differential equation for the master integral /, we use the 
differential equation (j7]) in variable t — u\ 



dg 1 da , 
a Hi + 2 ~dt 9 + Wl] U3 ' W3] W2 ' U2 ' = 



(11) 



and perform the inverse Laplace transform in t. This differential equation, using the following 
new parameters which are adapted to the symmetry of the problem 



W2 = w + x, ws = w — x, ii2 = u — y, W3 = u + y, 



takes the form 
where 

a 

(TO 



d 2 d 2 



d d 



dr z dr A dr dr 



a + t 2 CT 2 + t 4 (T4, 



(12) 



(13) 



(14) 



w 1 —2w l {u + w + x +y) + lQuwxy 



w\ + w\ 22 + CX20, 



(u + w — x — y) (u — w + x — y) (u — w — x + y) (u + w + x + y) 
+16 (w x — uy) {u x — w y) (uw — x y) 
w\oqi + coo, 



and F(r) = F Ul {r) is presented in Appendix A. The differential Eq. ( TT3l) is supplemented 
by the boundary conditions, namely f(r) vanishes at small and large r. 

Taking the occasion, we present here, for later use, a differential equation which is ob- 
tained from the inverse Laplace transform of Eq. (|T0|) 

(^-»^)^ + (^ + -^)^) + (- + ^^)/-(^-^«, d5) 
where F Wl (r) is given in Appendix A, and that for an arbitrary parameter a = u,w,x,y 

, 1 ih -> ^ - " 3 

2 da u da 



2^ + ^)'>> + ^ (1, M 



-FJr). 



(16) 



III. POISSON REPRESENTATION 



We derive here the integral representation for solutions of the homogeneous and the 
inhomogeneous differential equation flT3|) . which is the analog of the Poisson representation 
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for Bessel functions 18j. a is a quadratic polynomial in t 2 , so it has four zeros ±t a , ±tb- 
Assuming tf, > t a > 0, the four solutions of the homogeneous equation 

d d 



d 2 d 2 

°4 r ~n + °2 ~r r — + cr r 
dr z dr z dr dr 



f(r) = 0, 



arc 



-t b e tr 



: dt, 



-ta gtr 



dt, 



ta 



J. r 



: dt, 



h e tr 



Let us prove it, as an example, for the last term 



dt. 



(17) 



(18) 



d 2 



d d 



dr z dr 1 dr dr 



t,, 



t b 



dt 



dt 



d 2 



dt —= 

V a ( t ) 

d d 



dr 2 dr 2 



0-2 -r- r ~r + °o r 

dr dr 



t r 



ta ^/W) L 

i r tb dt 

2 



a 4 t —t + &2t — t + a — 
dt dt dt 



J r 



J r 



(19) 



•t a y/a(t) L 

where the last equation holds because a(t) = a + a 2 t 2 + a^t 4 . One integrates by parts, 
boundary terms vanish because a(t a ) = a(tb) = 0, and one obtains 



1 



1 



dte 



t r 



d 



d 



dt^ + ^dt^T) 



dte 



t r 



1 



1-1} 



0. 



(20) 



The solution of the inhomogeneous fundamental differential equation (|T3|) is obtained by 
noting that the integration contour in the inverse Laplace transform can be deformed to 
encircle all branch cuts on the left side of the complex plain 



f(r) 



1 [~°° 

— / dte tr \g(t + ie) - g(t-ie 
"ki Jo 



(21) 



The analytic properties of g have been analyzed by Fromm and Hill in 12] . The function g(t) 
has four overlapping branch cuts, similarly to the function P, in variable t on the negative 
real axis starting at 



h= u 3 + w 2 
-t 2 = u 2 + w 3 

-t 3 = U3 + WX + W3 
-t 4 = U 2 + Wi + w 2 



u + y + w + x, 
u — y + w — x, 
u + y + w — x + wi, 
u — y + w + x + Wi, 



(22) 



correspondingly. The fact that all U with i = 1,2,3,4 are not positive comes from the 
requirement that the integral in Eq. ([2]) is finite for positive values of Ui and W{. Fromm and 
Hill found the imaginary part of g for the particular ordering of t , . Their result is generalized 
here to an arbitrary ordering of £j and takes the form 

— [git + ie)-g(t-i e)] = — = [9{t x - t) In |/3 ,o| + B{t 2 - t) In |/3 3 , 3 | 

Z 7T Z Z \/ (7 



where 



-0(t 3 - t) In |/3 3 ,i &, 3 | - 0(t 4 - *) In \Po,i M] , (23) 



and Ccf. Ref. 



ftj = (24) 



12| with di = u 3 , a 2 = vo 2 , a 3 = t, a 23 = w 3 , a 31 = u 2 , a l2 = w x ) 



7o i0 = 2u 2 vo 1 vo 3 + (u 2 - u 3 + w x ) u> 3 + wi (—t + u 2 + vo 3 ) + u 2 (w{ - vo 2 + w 3 ), 

7 3)3 = 2u 3 vo 1 vo 2 + {-u\ + u\ + w\) w 2 + wi (-t 2 + ul + vol) + u 3 {w\ + wj- vol), 

7o,i = -2tw 2 vo 3 - (t 2 - u\ + vol) w s + w 2 (t 2 - uj + vol) + t (~ w l + W 2 + ^3)3 

73,1 = -2tvo 2 w 3 + (t 2 - ul + vol) w 3 — w 2 (t 2 — u l + w l) + t (-vol + w l + W D- (25) 

There is some arbitrariness in the form of Eq. (T231 as products of can be expressed in 
several ways 



Pfl,flPfl,U 0V,U fin, /J, ftp,P fip,Gl 

fifi,<j fiu,a ftp,® = 0cr,a-> (26) 

for {fi,u,p,a} being an arbitrary permutation of {0,1,2,3}. One can verify that Im(g) in 
Eq. (j23p . satisfies the differential equation ffTTj) with Im(P) as an inhomogeneous term. As 
a result, f(r) can be expressed in terms of the one- dimensional integral 

f(r) = / dte* r _=[0(ti-t) ln|^ ,o|+^2-t) ln|#, )3 | 
Jo 2 Vcr 

-0(t 3 - *) In |/? 3 ,i /?3, 3 | - 0(* 4 - *) m [/0o,i A,,„l] • (27) 

This is a main result obtained in this work. All properties of the master integral / can in 
principle be obtained from this integral representation. Some of them, however, can be more 
easily obtained from the differential equation ({TBI . The rest of the paper will be devoted 
to various properties and various ways of numerical calculation of f(r) in Eq. (T27]) and its 
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derivatives with respect to nonlinear parameters. This form of Eq. (|27|) is suited for the 
direct numerical quadrature, once we know all the singularities on the integration path, and 
this is presented in Appendix B. 

It would be interesting to investigate various expansions of f(r) in Eq. ( 127)) . in particular 
the expansion around w\ = is studied in Sec. VI. Here, we obtain the master integral / at 
w± = in a simple form which is convenient for the numerical evaluation. The a polynomial 
reads 



o"Ui=o = {u\ - u\ + w\ - w\) (u\ w\ - u\ wl) - t 2 (u\ - u\) [wl - wl) 
= 16 (w x — u y) [u x — w y) {u w — x y) + 16 1 2 u w x y, 



and the differential equation becomes [1 



16 u w x y — r - — h 16 r (w x — u y) (u x — w y) (u w — x y) 
dr dr 



where the inhomogeneous term F(r) using Eq.( lA5j) is 



f(r) = F(r) 



«>i=0' 



(28) 



(29) 



F( r )\ Wl - — 2 (uw x — uwy — uxy + w xy) Fi_ + 2(uw x — uwy + uxy — w xy) F 2 - 
—2{uwx + uwy + uxy + wxy) F 3 _ + 2 (uw x + uwy — uxy — w xy) F 4 _, 

(30) 



and where 

F2± 
F3± 

F,± 



Ei(-2 r u) e r ^~ w ~ x +y) ± Ei(-2 r w) e 



r (—u+w+x—y) 



Ei(-2rw) e »- (-«+»-*+») ± Eif-2 



rule 



r (u—w+x—y) 



Ei(2 rx) + Ei(2 ry) - Ei(2 r (x + y)) + In 
±Ei(-2 r (u + w)) e T ( u+w+x+ y\ 



Ei(-2 r x) + Ei(-2 ry) - Ei(-2 r (x + y)) + In 



uw (x + y) 
xy (u + w) 

uw (x + y) 



-r (u+w+x+y) 



xy (u + w) 



-r (u+w—x—y) 



±Ei(-2r(u + tw))e 



r (u+w—x—y) 



(31) 



with Ei being the exponential integral function. The solution of this differential equation 
assuming p 2 > 0) in terms of solutions of the homogeneous equation was presented in Ref. 



f(r) 



1 

IQuwxy 



roc pr 

I (pr) / dr' F(r')K (pr') + K (pr) / dr' F(r) I (pr') 

Jr JO 



(32) 



where 



V 



(w x — uy) (wy — u x) (u w — x y) 



(33) 



uwxy 

and Iq, K are modified Bessel functions. A more convenient, however is the integral form 
ofEq. (Ej 



f(r) 



— oo r — oo 

ln|/3 ,o|+ / ln|/3 3 , 3 | 

u—w—x—y J — u — w-\-x-\-y 

— oo r — oo 

hi|/3 3 ,i/3 3 ,3|- / ln|/3 ,i/?o,o| 

u—w+x—y J —u—w—x+y 



J, r 



= dt, 

2v^ 



(34) 



with the following 7^ coefficients 

7 0j o = —4:(uwx + uwy — uxy — wxy), 

73,3 = 4:(uwx + uwy + uxy + wxy), 

r(/3 3i i /3 3j3 ) = — 16uwxy/(t + u + w — x + y) + 4 (uwx — 111111/ + uxy + wxy), 

r(A),iA),o) = —16uwxy/(t + u + w + x — y) — 4:(uwx — uwy — uxy — wxy), (35) 

where T is defined in Eq. (|B3]). In the particular case of exchange integral w = u, x = y, 
o"oo vanishes and the master integral becomes 



f(r) 



■ w x 



dt In 



2 (w+x) 



t + 2(w-x) 



t-2(w-x) 



dt In 



2 (w-x) 



t-2(w + x) 



dt In 



2 w 



t + 2(w + x 
t + 2(w-x 



A r 



t + 2(w + x) 



J r 



(36) 



An equivalent integral representation, obtained from the differential equation (|29|) was pre- 



sented in Ref. 



13| . but this form is more convenient for the numerical evaluation. 



IV. TAYLOR EXPANSION IN r 



Here, we find the Taylor expansion in the internuclear distance r of the master integral 
/. As it was noticed by Fromm and Hill in [l^ this expansion is absolutely convergent, 
therefore it is another way to calculate the master integral as well as its derivatives with 
respect to it, and W{. Following Ref. 12|, we find the initial terms and the recurrence 
relation for subsequent terms of the Taylor series. 



The expansion in small r of the master integral has the form 

oo 

f(r) = J2r k [fi 1] ln(r) + /f]. 



(37) 



fc=i 



The leading term, can be obtained by noting that the limit r — > corresponds to the 
one-center integral 12 1 



1 



16 7T 2 



e -ori-/3r2-7ri2 

<fVi / d 3 r 2 5-5 



1 

: 2^ 



6 2 1/3 + 7 



+ L 2 



a + 
a + 7 



L 2 



/3 + a 
^+7 



(38) 



where L2 is the dilogarthmic function [181 ] . Let us introduce for the later use the following 
Xi symbols 



X 

x 2 



2wi 



7l~ 



— + - In 2 

6 2 



2 u + u>i 
2w + wi 



+ L 2 1- 



2 (u + w) 
2u + w\ 



L; 



2(u + w) 
2 w + w\ 



In 



In 



2 u + iox 
2 (u + iu) 

2 (it + w) 



X r = - ln[r 2 (2M + w 1 )(2w + Wi)] + 7 E . 



(39) 



where 7e is the Euler constant. The leading terms in small r expansion using Eq. (1381) are 



0, 

x . 



(40) 



and /o = f_{ = 0. The next terms of this expansion f 2 , f 2 2 ^ are obtained from the large 



t asymptotics of the function g 



^ = 72" + 73 + °( t )■ 



t 2 ' t 3 

Since g(t) satisfies the differential equation ( TTTT) with the inhomogenous term P(t) 

P(t,Wi,U3,W 3 ,W 2 ,U 2 ) 



(41) 



P(t) = w\ [2 + ln(u 2 + u 3 + wt) + ln(wi + w 2 + w 3 ) - 2 \n(t)] + o^ 1 )., 



(42) 



the 1/t 3 coefficients are 



-2, 



7(2) 



ln(w 2 + u 3 + Wi) + ln(wi + w 2 + W3). 



(43) 
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As noticed by Fromm and Hill [12| . they are related to coefficients by 

-a) _ _^ 

h - 2 , 



(44) 



where ip is the Euler -^-function 18] and i[)(3) = — 7e + 1 + 1/2. The next terms in r expansion 
can be obtained from the fundamental differential equation ( fl~3l) . As the inhomogeneous term 
F(r) has a similar expansion 



F(r)= 5> fc (if> Hr) + Fi*>): 



(45) 



k=-l 



the recurrence relations are the following 



f (i) 



(2) 
fc+1 



a 4 (Jfe + 1) k 2 (k - 1 



A-.-^o/ii ; 3-^(*-l) a /iI ; 1 -^2(*-l)/ J SL ) 1 -a42fc[2(fc + l) (*-!) + !] 



(46) 
(i) 



a 4 (A; + 1) k 2 (k - 1) 

The resulting first terms in the expansion of the master integral / in powers of r are 

3* 



/(r) = rX + r 2 (x r --) 



+r" 



w\ X wi (u 2 + w 2 + x 2 + y 2 ) X - w Xx - u X 2 - 2 (u + w) 

12 4~ + 6 
it/ 2i?/(2mwX + hX 1 +wX 2 ) 



3 u>i 



54 ^ 



wf + 



3w 2 

W\{u + w) 3uw — 3xy — 11 u 2 — 11 w 2 — 9x 2 — 9y 2 



12 



+ 



36 



wf u z + w z + ar + y z 

+A r I H 

1 18 6 



0(r 5 ). 



(47) 



Let's consider now the special case of W\ = 0. The recurrence relation takes the form: 

_ f (i) 



f (l) 

Jk+l ~ a 20 (k + l) 2 



f (2) 
Jk+l 



-Pfc 2) - Opo/fc^i _ 2 (i) 

a 20 (fc + l) 2 fc+l /fc+1 ' 



(2) 



(48) 
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and the first terms in r-expansion are 

1 



f(r)\ 



toi=0 



1 i u \ 1 i (-, w 

— In 1 + - +— In 1 + - 

2u V w/ 2w V u 



+ r 



3 1 
— - + 7 + - ln(4 r 2 u w) 



u w xy xy 
3 ~ ¥ + 18m + 18w 



m w 2 + x 2 + y 2 w xy 
4 + 12u 18 m 2 



ln(l + - 

w 



w u 2 + x 2 + y 2 uxy 
1 + 12tu 18 



In 1 



+r" 



uw — xy w 2 + u 2 l,o 9 9 9s / 3 1 , / , 9 

^— + g (n 2 + w 2 + x 2 + y 2 ) ( -- + 7+2 ln(4r 2 uw; 



12 



+0(r 5 ). 



(49) 



In order to obtain the Taylor expansion of f(r, n) 

oo 

f(r,n) = J2r k ln(r) + if (n)] . 



(50) 



fc=i 



one expands Eq. (fl~5l) in power series in r and obtains the following recursions for coefficients 



(i) 

k+2 



n 



aoo fi 1] (n) + cr 02 (n - 1) »f£>(n - 2) 



2 /(I)/ 



<r 20 (A; + 2) (A; + 1) 
+( x 22 (n-l) 2 (A; + 2) (* + 1) /« 2 (n - 2) 

+(n - l) 2 (A; + 4) (k + 3) (A; + 2) (fc + 1) fj^(n - 2) 

+(„ _ 3 ) (n - 2) 2 (n — 1) (A; + 2) (A; + 1) /« (n - 4) 



(51) 



(2) 



(T 2 o (k + 2)(k+ 1) 

,(2) 



i^LfcCn - 1) - *20 (2k + 3)1 



WrW-^o 2 (n-l) 2 /f (n-2) 



(i) 

fc+2 



n 



a 22 (n - l) 2 (* + 2) (fc + 1) fi%(n - 2) + (2 fc + 3) /# 2 (n - 2) 



(i) 



_(„ _ 3 ) ( n - 2) 2 (n - 1) [(* + 2) (k + 1) /g 2 (n - 4) + (2 fc + 3) - I ) 

-(n - l) 2 [(* + 4) (fc + 3) (A; + 2)(k + 1) /i 2) 4 (n - 2) 



+2 (2 fc + 5) (A; 2 + 5 k + 5) /f 4 (n - 2) 



(52) 



where (0) = . The usefulness of the calculation of /(r) and f(r,n) via Taylor series 
in r and the above recursions needs to be verified numerically. Nevertheless, the expansion 
terms say a lot about analytic properties of /, as a function of Wi, x, y, u and w. In particular, 
the expansion in Eq. ( 1471) suggests that / is an entire function in x,y at fixed r, u, w, w\. 
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V. GENERAL RECURSION RELATIONS 



In order to use the exponentially correlated basis in molecular calculations, one has to be 
able to calculate derivatives of the master integral / with respect to nonlinear parameters. 
Differential equations (I13|ll5p and (|TB|) are used below, to express the derivative of /(r) with 
respect to an arbitrary parameter a = wi, x, y, u, w in terms of f(r), f'(r), f"(r), and f"'(r). 
This is done as follows. From Eq. ffl3l) . the fourth derivative f^(r) can be expressed in 
terms of lower derivatives of f(r) and we use this in Eqs. (fl5|) and (Tl6|) to eliminate fourthth 
derivative. Next, differentiate these equations with respect to r and again eliminate fourth 
derivative. Doing this differentiation and elimination three times, one obtains four equations 
for four unknowns: the derivative of / with respect to Wi, and its three further derivatives 
with respect to r. The solution of these linear equations for a = w\ is 

df(r) _ <9cx n , / _ dap _ da? 



dwi 



2 5 a =r6^f\r)+w 1 (-2a a 2 -w 1 a 2 ^ + 2w 1 a ^) [f"(r) + r f"{r)] 



da 
dwi 



dw\ 



a ( 4 u>i a + 2 w\ 



da 
dw\ 



0"2 



da 2 
dwi 



[f(r)+rf(r)] + ... 



(53) 



and for a = u,w,x,y 
df(r) 



28 a 



da 



^ ^\lHr)+rf(r)] + ..., 



ap[2w l 



da da 

where . . . denotes inhomogeneous terms presented in Eqs. (IC3IIC5I) . and where 



(54) 



5 = 4 w\ a — a ^ 



Au 1 



Aw 1 



Ax 1 



(55) 



The solution for higher order derivatives with respect to r can be represented in the matrix 
form. Let us introduce a symbol / to denote 



f(r) 



( f(r)\ 
f(r) 
f"(r) 
\ f"'(r) J 



(56) 



Then 



da <*> 



(57) 
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where A a and F a are presented in Appendix C. These vector differential equations with 
respect to the nonlinear parameter a can be used as recursion relations to express the 
integral with any powers of electron distances in terms of the master integral f(r) and it's 
first three derivatives with respect to r. These recursions become unstable when w\ ~ 0, 
5 ~ or o"o ~ 0. These cases require a separate analysis. In the next Section we consider 
recursion relations for W\ = which has applications for standard molecular integrals. 

The particular case of <j$ = happens for an exchange integral where w = u and x = y. 
Here, one can use recurrence only for f'(r) to generate the analytic expressions for derivatives 
with respect to w — u and x — y. They are quite compact as there will be only three 
independent parameters, and can be expressed in terms of f(r) and f'(r). The derivatives 
with respect to remaining nonlinear parameters can be calculated afterwards, by setting 
w = u, x = y. 



VI. RECURSIONS FOR STANDARD MOLECULAR INTEGRALS 



We describe in this Section the evaluation of standard molecular integrals, where w\ = 0. 



Recursion relations for the integral f{r,n) in Eq. ([3]) have been presented in Ref. 13 ]. 
Here we rederive them using a more compact notation. In particular cases, close to spuri- 
ous singularities, we propose to use these recursions to derive an analytical expression for 
individual integrals, which are then implemented in the numerical codes. 

Let us take Eqs. ( fl3l) and ( Ti~5l) . and differentiate them with respect to W\, n and n — 1 
times respectively at wi = 

r ooo f(r, n) + a 20 f'(r, n) + r cr 20 f"{r, n) + (n - 1) n r <t 02 f(r, n - 2) 
+ (n - 1) n a 22 f'(r, n - 2) + (n - 1) n r o 22 f"(r, n - 2) + 2 (n - 1) n f {3) (r, n - 2) 
+ (n - 1) nr / (4) (r, n - 2) + (n - 3) (n - 2) (n - 1) n f'(r, n - 4) 
+ (n-3)(7i-2)(n-l)nr/"(r,n-4) = F Ul (r,n), (58) 

<roo f(r, n) + a 20 f"(r, n) + (n- l) 2 a 02 f(r, n - 2) + (n - l) 2 a 22 f"(r, n - 2) 
+ (n - l) 2 / (4) (r, n - 2) + (n - 3) (n - 2) 2 (n - 1) /"(r, n - 4) = F Wl (r, n - 1) . (59) 

These are two linear equations for three unknowns f(r, n), f'(r, n), f"(r, n). The third equa- 
tion is obtained by elimination of f"(r,n) and further differentiation with respect to r. The 
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solution of these three equations for f(r,n) and f'(r,n) are 
f(r,n) 



[_(„ _ 2 ) {n - 1) a 02 f(r, n - 2) + (n - 1) r a 02 f(r, n - 2) 



f(r,n) 



o"oo 

+2 (n - 1) a 22 /"(r, n - 2) + (n - 1) r a 22 / (3) (r, n - 2) 
+(n - 1) (n + 2) / (4) (r, n - 2) + (n - 1) r / (5) (r, n - 2) 

+4 (n - 3) (n - 2) (n - 1) f'(r, n - 4) + 2 (n - 3) (n - 2) (n - 1) r / (3) (r, n - 4) 

+2 F W1 (r, n - 1) + r F' Wi (r, n - 1) - < (r, n)] , (60) 

— — [(n - 1) r (T 02 /(r, n - 2) + (n - 3) In - 2) (n - 1) n /'(r, n - 4) 

+(n - 1) n a 22 /'(r, n - 2) + 2 (n - 3) (n - 2) (n - 1) r /"(r, n - 4) 

+ (n- l)ra 22 /"(r, n - 2) + 2 (n - 1) n / (3) (r, n - 2) + (n - 1) r / (4) (r, n - 2) 

+rF ttl (r,n-l)-F ul (r,7i)], (61) 



where 



F x (r, n) 



d" 



F x (r) 



(62) 



for X = Wi, U\. Equation ( |60l) allows one to obtain integral f(r, n) with an arbitrary power 
n — 1 of ri2 in terms of /(r), for example 



f(r,0) 
f(r,l) 
f(r,2) 



f(r) 



4 



h (ru) h (rw)j (r x)j {ry), 



— [r ao2 f'(r) + 2 a 22 f"(r) + r a 22 f^\r) + 4 /^(r) + r /( 5 )(r) 



+2F Wl (r,l)+ri^ 1 (r,l)-<(r,2)]. 



/(r,3) = 



3h 1 (ru)h 1 (rw)j 1 (rx)j 1 (ry)-h (ru)h (rw) j (r x) j (r y) 
+h (rw) h 2 (r u) j {r x) j (r y) + h (ru) h 2 {r w) j (r x) j (r y) 
+h (ru) h (r w) j 2 (r x) j (r y) + h (ru) h (r w) j (r x) j 2 {r y) 

where j n and h n are modified spherical Bessel functions 

1 d \ n sinh(a;) 



(63) 



Jn{X) 

h n (x) 



X 



X' 



1 d \ n exp(— x) 
x dx, 



x 



(64) 



One notices, that f{r,n) for odd n can be expressed in terms of j n and h n only, and their 
numerical evaluation is straightforward. In contrast, numerical evaluation of the formulas 
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for even n is much more difficult. If 020 is not small, then derivatives of f(r) can be reduced 
to /(r) and f'(r) using 

/W(r) = - J_ 
If (J20 is small, then the inverse recursion 



<7 Q0 [(n-2) / (n " 3) (r) + r / ( - 2) (r)]+a 20 (n-1) ^(r)- (r)] (65) 



/W(r) = - ? [r a 00 / (n+1) (r) + (n + 2) ffao /( n+2 )(r) + r a 20 / {n+3) (r) - F% +1 \r)] 

(n + 1) (Too 

(66) 

is stable and can be used to obtain all the derivatives including the function /(r) itself. 
Alternatively f^ n \r) can be obtained directly, using the integral form of Eq. ( )34|) . If o"oo is 
small but not 020, then these formulas can be rewritten using Eq. (IBTj) to the form without 
(Too in the denominator, for example 

f(r,2) = — [a 02 f(- 2 \r)-ra 02 f- 1 \r)-a 22 f(r)-rcr 22 f\r)-3r(r)-rr(r)}+... (67) 
020 

where f^~ n \r) is defined by 

/(-™)(r) = / f(- n+1 \r)dr, (68) 







and can be obtained from Eqs. ( J34l) or fl66l) . It may lead however to other type of numerical 
instabilities, if the branch point starts at 0, see Eq. (|22p . and in this case one can build 



the explicit table of integrals, as in Ref. [14], ll5[. When both <j 2 q ~ and cr 00 w 0, then 



parameters x, y are small and in this case the Neumann expansion [5| can be applied. 

What remains are the derivative of / with respect to a = u,w,x,y at W\ = 0. We 
adapt here the derivation of corresponding formulas from Ref. [ljj. One takes fundamental 
differential equations in u\ and a at w\ = 

r ^00 f(r) + 2 o f\r) + r a 20 f"(r) = F Ul (r) (69) 

1 da 00 df(r) 1 da 20 df"(r) 

2 &r /(r) + a °° -^r + 2 (r) + " 20 -*r = ~ Fa{r) (70) 

differentiates the first equation with respect to a, eliminates df"(r)/da and differentiate 
resulting equation again with respect to r. The obtained equation for the derivative of /, 
using da 20 /da = <7 20 /a is 

df I r(l da 00 l\ G a (r) 

77T = -77- + o — -KT--,}f( r ) + — - > ( 71 ) 



9a 2 a 2 \ <t o 9a a / <j o 
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where 



and 



GJr) 



1 

2a 







d 



da )K(r)-(2 + rf r \K,n 



- w\=0 



G w {r) = — (uy - wx) (ux - wy) (-F 1+ - F 2+ + F 3+ + F 4+ ) 



+(u + w) (x — y) 
+(u — w) (x + y) 



GJr) 



1 

u 



(uy - wx) (ux - wy) (-F 1+ - F 2+ + F 3+ + F A+ 



+ (u + w) (x — y) 
+ (u — w) (x + y) 



GJr) 



1 

x 



(uy — wx) (ux — wy) (F] 



+ (u + w) (x — y) 
+ (u — w) (x + y) 



G y (r) 



1 r 

y 



(uy — wx) (ux — wy) (F] 



+ (u + w) (x — y) 
+ (u — w) (x + y) 



uw — xy) 
uw — xy) 



uw — xy) 
uw — xy) 



uw — xy) 
uw — xy) 



uw — xy) 
uw — xy) 



—Fi + + F 2H 
-^ 3+ + F iA 



F\+ — F 2+ ) 
-F 3+ + F i+ 

f + F 2+ — F 3+ 
-F 3+ + F 4+ ) 
—Fi + + F 2+ ) 

\- + F 2+ — F 3+ 
F 3 + — F i+ ) 
F\+ — F 2+ 



(72) 



(73) 



Here <Jq from the denominator can be removed, when necessary, as in the case of f(r,n) by 
introducing f^~ n \r). The formula ( 17TI) allows one to obtain integrals with arbitrary pow- 
ers of electronic distances. It requires however, further investigations of various numerical 
instabilities when the denominator approaches 0. In particular, a Q0 vanishes for exchange 
integrals. In this case w = u,x = y and one can generate a table of integrals avoiding 
recursions, by Taylor expansion of the fundamental solution in Eq. (I2~T]) in w±, w — u, and 
x — y. Alternatively, one can take recursion relations presented in Eq. f l6T|) and integrate it 
once over r. 



VII. SUMMARY 

The evaluation of two-center two-electron integrals is a difficult task. Depending on a 
physical problem different approaches are employed. For the calculation of nonrelativistic 
energies of a two-electron diatomic molecule the best way is by the use of Kolos-Wolniewicz 
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basis functions, namely w\ = and arbitrary nonlinear parameters u,w,x,y. The master 
integral can be calculated numerically using Eqs. ( |34|) and ( 157)) . Integrals with powers of 
interparticle distances can be obtained by using analytic recursions Eq. (16TJ)) . ( IBS)) , and ITT)) 
followed by numerical evaluation of obtained formulas. 

For the calculation of QED effects in a two-electron diatomic molecule, the best way is 
to use the most general basis functions with nonvanishing w±, u, w, x, y. The master integral 
can be evaluated according to Eq. (1271) . and integrals with additional positive powers of 
interparticle distance can be obtained through recursion relations Eq. (1571) . Integrals with 
negative powers can be obtained by numerical integration with respect to the corresponding 
nonlinear parameter. In order to demonstrate correctness of obtained formulas for two- 
center two-electron integrals, Table |TJ presents a simple, examplary calculations of the ground 
electronic state of the H2 molecule. 

TABLE I: Examplary calculations of the Born-Oppenheimer energy at the equilibrium distance 
R = 1.4 of the ground state of H2 molecule, obtained using functions of the form Eq. (JT]), with 
parameters w\ = — 0.5, u = 1, w = l,y = 0.125, x = —0.125 defined in Eqs. f)12[) . TV is a number of 
basis functions, and n i + ni + n.3 + 714 + 715 < Q. The result with 00 is accurate to all digits and 
was obtained in 
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14] 



n 


N 


.EboO-'"-) 





1 


-1.038 401456 


1 


4 


-1.157108 377 


2 


13 


-1.172 368483 


3 


32 


-1.174 254201 


00 




-1.174475 714 



An additional approach to evaluate two-center integrals is by the Taylor expansion in 
small internuclear distance r. Its applicability has to be verified numerically. Nevertheless, 
it serves at present as a simple check of various formulas for recursion relations and the 
master integral. The advantage of calculations by Taylor series, is that expansion terms are 
quite simple. They are related to the known one-center (helium-like) integrals and can be 
obtained recursively ( 146) I5T) 152)) . A similar Taylor expansion can in principle be obtained 
for many-body three-linked integrals (at most three odd powers of interparticle integrals). 
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This would open a window for high precision results in many electron diatomic molecules. 
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Appendix A: Inhomogeneous terms 



The first order differential equation ([8]), satisfied by the general four body integral, defined 
by Eq. (3), involves the inhomogeneous term P a for a = u^Wi. They are related to each 
other by 

P Wl = P(w 1 ,ui]w 2 ,u 2 ;w 3 ,u 3 ) 

= P(wi,Ui]W 3 ,U 3 ]W 2 ,U 2 ) 
P U1 = P(u 1 ,W 1 ]W 2 ,U 2 ]U 3 ,W 3 ) 

P(w 2 ,u 2 ;w 3 ,u 3 ;w 1 ,u 1 ) 
P(u 2 ,w 2 ;w 3 ,u 3 ;u 1 ,w 1 ) 
P(w 3 ,u 3 ;w 1 ,u 1 ;w 2 ,u 2 ) 

P(u 3 ,w 3 ;u 1 ,w 1 ;w 2 ,u 2 ) (Al) 



11' 2 



a 2 



W3 



The expression for P was obtained in [13J and is the following 



In 



P(lUi,U 1 ]W 2 ,U 2 ]W 3 ,U 3 ) 

UiW l [(ui + w 2 ) 2 - uj] 

(-«! + U 3 - W 2 ) {Ui +U 3 + W 2 ) 

| u 1 w 1 [(u 1 + u 3 ) 2 -wj] 

(-Mi - U 3 + W 2 ) (Ui + U 3 + W 2 ) 



111 



u 2 + U 3 + Wi 
U\ + U 2 + Wi + W 2 

Wi+W 2 + w 3 



ui w\ + u\w\ — u\ w 3 + Wi w 2 {u\ + u\ — w 



Ul + U 3 + Wi + w 3 _ 

Ui +u 2 + w 3 



111 



111 



{—Wi -w 2 + W 3 ) (Wi + w 2 + w 3 ) 

u\ w\ — ulw 2 + u\ w| + Wi w 3 (uf + u 2 3 — w 2 ) 
{—Wi + w 2 - W 3 ) (Wi + w 2 + w 3 ) 

u 2 (u 2 + wx) {u\ + uj- wl) - u\ {u{ + uj- wf) 
{-u 2 + u 3 - vox) (u 2 + U 3 + Wx) 

u 3 (u 3 + wx) (u\ + ul~ wl) - uj (uj + uj- wl) 
(u 2 -u 3 - Wi) (u 2 + U 3 + Wi) 

wi [w 2 (uf -ul + wl) + w 3 (ul -uj + wl)} 



111 



In 



Ul + U 2 + Wi + w 2 

Ui+U 3 + w 2 
Ul + U 3 + Wi + w 3 

Ui+U 3 + w 2 



Ul + U 2 + Wi + w 2 
Ui+U 2 + w 3 



(wi —w 2 — w 3 ) (wi +w 2 + w 3 ) 

wi [u 2 (ul + u 2 3 - wl) + u 3 (ul + w\- wl)] 
(-u 2 -u 3 + Wi) (u 2 + U 3 + Wi) 



111 



111 



Ui +U 3 + Wi 
U 2 + U 3 + Wi 



w 3 



u 2 + u 3 + w 2 + w 3 

Wi+W 2 + w 3 

u 2 + u 3 + w 2 + w 3 



(A2) 
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In the case the nonlinear parameter a is the combination of Wi,Ui as in Eq. ( fl2|) . then P a 
is given by 



p 

p x 
p 



p 4- p 

P — P 

1 W2 1 IDs 

P + P 



P„ 



P — P 



(A3) 



The inhomogeneous term in the fundamental differential equation (1131) is the inverse Laplace 
transform of P 



1 



2ni 



dte tr P, 



U\=t 



(A4) 



In particular, F Ul is given by 
F Ul (r) = 



Wl{ - + 2 w 1 + u + w-x + y ^ e _ r{u+w+Wl _ x+y) 
_„,, , 1 + u + w + x + y ^ ( _,.,„+„.+,._,„ 



+ 



x 



+ 



X 



w 



- (u— w — x + y) + 2uw(x — y) + 2xy (w — u) 



Ei(-r (2U + Wl )) e r(u- W -x+y) _ Ei (_ r ( 2w + w ^ e -r («- 



— 3-( , u + w + x + y) + 2-uw(x + ?/) + 2x?/('U + u7) 
(w{2 r (x + y)) - Ei(-r (wi - 2 x)) - Ei(-r (w x - 2 y)) 



In 



(x + y) (w\ + 2u) (wi + 2w) 
_(u + w) (wi — 2 x) (wi — 2 y) 

+ < x — > —x, y — y 



e -r(u+u,+x+y) + Ei (_ 2r (U + V))) e T {u+W+X+v) 

(A5) 
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where Ei is the exponential integral function, and F Wl is 



r {u + w + x + y) + 1 



Wi 



Wi 



Wi 



Wi 



2u + wi 2w + wi w\ — 2x w\ — 2y 



^2 f {u + w) (u + x) (u + y) (u + w) (w + x) (w + y) 
r \ 2u + Wi 2w + wi 

(u + x) (w + x) {x + y) (u + y) (w + y) (x + y) {u + w + x + y)" 



Wi — 2 x 
r (u + w + Wi — x + y) + 1 



w x -2y 



+ 



Wi 



+ 



Wi 



+ 



Wi 



+ 



Wi 



2u + wi 2w + Wi w x — 2x Wi + 2y 

2 



Wl 1 . , 

+ — H — w l + A(xy + ux + wx — uy — wy — uw) + (wi + u + w — x + y) 



2 {u — w) {u + x) {u — y) 2{u — w) {w + x) {w — y) 



2 u + Wi 



2 w + Wi 



2 (u + x) (w + x) (x + y) ^2 (u — y)(w — y) {—x — y) 



wi — 2 x 
+ \ x — > —x, y — > —y 



Wi + 2y 



-r (u+w+wi—x+y) 



(A6) 



One notes, that F W1 involves only the exponential and rational functions. 



Appendix B: Numerical evaluation of the master integral 

For the numerical evaluation of the master integral in Eq. ff27j) . one should find all 
singularities on the integration path 12J. For this one decomposes a — 7? - into products 

a - 7q = (u 2 -u 3 + wi) (u 2 + u 3 + wi) (t-u 2 - w 3 ) 

(t + u 2 + w 3 ) (wi - w 2 + w 3 ) (wi + w 2 + W3) 
a - 7J3 = (u 2 -u 3 - wi) (u 2 + u 3 + wi) (-t + u 3 + w 2 ) 

(t + u 3 + w 2 ) (wi + w 2 - W 3 ) (Wi +w 2 + w 3 ) 
a- - 7^ = (t + u 3 - w 2 ) (-t + u 3 + w 2 ) (-t + u 2 - w 3 ) 

(wi +w 2 - w 3 ) (t + u 2 + w 3 ) (wi -w 2 + w 3 ) 
- 7o,i = (~t + u 3 - w 2 ) (t + u 3 + w 2 ) (t + u 2 - w 3 ) 

(wi +w 2 - w 3 ) (—t + u 2 + w 3 ) (wi -w 2 + w 3 ) (Bl) 

and observes that /?o,o vanishes on the negative axis at t = t 2 , f3 3t3 at t = ti, f3 3t i at t = t 2 , 
and /3q,i at t — t\. All singularities at other points t — Ui—Wj, t + Ui — Wj and t — Ui + Wj cancel 
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out among various In /3's, because the original integral Eq. fl2]) is finite at these points and 
g is analytic. The cancellation can be achieved explicitly, by combining the corresponding 
/3's under the common logarithm. For this one has to consider all possible orderings of t^. 
Because of the symmetries lf)2 and A -h- B one can assume, without loosing generality, 
that ti > t 2 and t 3 > t 4 . We remain with six possible orderings of ti, t% > t 2 > t 3 > £4, t\ > 

t 3 >h> u, h>t 3 >u> t 2 , t 3 > h > t 2 > u, h>u>u> t 2 , t 3 >u>h> t 2 , 

and six corresponding representations of the master integral / 

(rt'z pt3 rti r— oo \ tr 

jf \n\/3 ,o\+J t ln|/3o,o/3 3 ,3| + ^ In |/?o,oAM " j[ In \/3 0>1 /3 3>1 \j -^dt 

In |/3 ,o| - / In \/3 0>2 \ + / In |A,, /Aj,i| - / In |/3 ,i /? 3 ,i| — ^ 
ti h-i Jt 2 Jt 4 J 



*3 ft4 fti /■— oo \ gtr 



ti Jt3 Jti Jt 2 



In |/3 ,o| — / In ^0,2! - / In |A>,oA>,i Ad ~ / hi |/? 0) i Aj.iI H= d< 



/>t2 rfc* /•— 00 \ tr 

/ ln|/3 ,o/3o l2 |- / ln|/3 ,2|+/ In |/W&,il ~ / h|A),iAj,il -pdt 



*3 •'ti •'ta J ti 

rt\ rt4 rt 2 r—oo \ t r 

/ In |^o,o /0o,2|- / ln|/3 0)2 |- / In |/3 , A),i A) >2 | - / In |/3 ,i /3 3 ,i| -7= dt 

Jt 3 Jtx Jt 4 Jt 2 / \° 

ti rti rt 2 p— oo \ 

hi 1/30,0/30,21- / ln|/3 2 /3 ,i/3o, 2 |- / ln|A) j0 /V A), 2 | - / ln|A>,i 

i 3 "><4 -><2 / 



In order to eliminate singularities algebraically, the products and ratios of can be combined 
together [t3] 

Vo- + 7 7' + y/a y/a + T 

where 

r = £±2?: ( B4) 

7 + 7 
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As a result, no singularity appears on the integration path, which is demonstrated by the 
following decompositions 

a - r 2 (/3 ,i = ~(t-u 3 - w 2 ) (t + u 3 - w 2 ) (t-u 3 + w 2 ) (t + u 3 + w 2 ) 

(t-u 2 - w 3 ) (t + u 2 - w 3 ) (t-u 2 + w 3 ) (t + u 2 + w 3 )/ (4 1 2 ) 
a - r 2 (/3 0i0/ %,i) = ~{u 2 -u 3 + Wt) (u 2 + u 3 + wi) (t-u 3 - w 2 ) (t + u 3 - w 2 ) (t - u 2 - w 3 ) 

(wi + w 2 - w 3 ) (t — u 2 + W 3 ) (tOi + w 2 + w 3 )/ (t — u 2 — Wi — w 2 ) 2 
a - r 2 (/3 0i0 ^0,2) = -(-«2 + U3 + wx) {u 2 + u 3 + wt) (t + u 3 - w 2 ) (t + u 3 + w 2 ) {t-u 2 + w 3 ) 

(t + u 2 + w 3 ) (wi -w 2 + w 3 ) (wi + w 2 + w 3 )/(t + U 3 + Wi + w 3 ) 2 
cr - r 2 (/3 ,o /3o,i) = ~{U2 -u 3 + IWi) (u 2 + U 3 + Wt) (t-u 3 + w 2 ) (t + u 3 + w 2 ) (t + u 2 - w 3 ) 

(w 1 + w 2 - w 3 ) (t + u 2 + w 3 ) (w l + w 2 + w 3 )/(t + u 2 + wi + w 2 ) 2 (B5) 

Table [TT] presents numerical results for / for several internuclear distances r and nonlinear 
parameters wi,x,y,u and w. For r < 0.01 / is evaluated from the Taylor series up to r 8 , 
and for r > 0.01 from the integral representation in Eq. (1B2|) . 

TABLE II: Master integral / for selected values of internuclear distances r and nonlinear parameters 
wx u 2 w 2 u 3 w 3 /(0.1) 10 2 /(1.0) 10 3 /(10) 10 13 

2.5 2.0 1.5 1.0 0.5 1.539 488 720 658182 3.811561883 331994 2.916 697700 943 504 
2.0 2.5 1.5 1.0 0.5 1.575 469 059 882 717 3.538 642196 033 083 1.056 827512 869 080 
1.5 2.0 2.5 1.0 0.5 1.592 898118 308 067 3.359 218 711378 032 0.983486 526 526 378 
1.0 2.0 1.5 2.5 0.5 1.687851128 764463 3.142 086 593 091502 0.648185 277118 373 
0.5 2.0 1.5 1.0 2.5 1.687626 825 828 684 2.985 309 587884137 0.621477317148 430 
-0.5 2.0 1.5 1.0 2.5 2.285 510 252 707 772 5.843 903 698 492 676 2.026400131640 827 
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Appendix C: Exact form of general recursion relations 



The recursion relations for the general master integral, as given by Eq. (|57|) involves the 
4 by 4 matrix A a . For a = x,y,u,w, A is the following 

/ _1 <M _r f dS_ _ 2_5 dao\ w\a a rwja a ^ 

da 



4 9a' 
2 ' 

2 ' 
r (7Q 95 

4 w'f da ' 



25 dap 
ao da 



raa 
2 
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2 ' 
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2 CT(J 



r 95 
4 9a 
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where 
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and the inhomogeneous term Fq, is 
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For a = Wi matrix A is 
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and the inhomogeneous term 



4 (to (^dwi ^1 dwi^j (^^1 ctq) drdw\ ao dr 3 ^1 CT0 dr^dwi 

V w\ J dr ^ drdwi 1 9r 3 1 dr^dwi 
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